La AMO è una oscillazione delle temperature superficiali nel Nord Atlantico, in particolare indica le oscillazioni ogni 60-80 anni che occorrono nel tratto di oceano compreso tra l'Equatore e la Groenlandia.
La regione della Groenlandia, soprattutto nell'ultimo decennio, è nota per essere molto influenzata dai cambiamenti climatici, l'aumento delle temperature e in particolare lo scioglimento dei ghiacciai può avere conseguenze devastanti per l'intero globo.
In questo notebook si vuole andare a studiare, se esiste una relazione tra temperature della Groenlandia e l'AMO index.
# importazione delle librerie
import os
import numpy as np
import pandas as pd
import geopandas as gpd
import datetime as dtm
from datetime import datetime
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import scipy as sp
from scipy import stats
import seaborn as sns
import xarray as xr
from PIL import Image
flist = [os.path.join(path, name) for path, subdirs, files in os.walk("./CRUTEM.4.6.0.0.station_files/") for name in files]
flist=flist[1:]
print ('\n',flist[0:5])
nst=len(flist)
['./CRUTEM.4.6.0.0.station_files/CRUTEM.4.6.0.0.station_files\\01\\010010', './CRUTEM.4.6.0.0.station_files/CRUTEM.4.6.0.0.station_files\\01\\010030', './CRUTEM.4.6.0.0.station_files/CRUTEM.4.6.0.0.station_files\\01\\010050', './CRUTEM.4.6.0.0.station_files/CRUTEM.4.6.0.0.station_files\\01\\010070', './CRUTEM.4.6.0.0.station_files/CRUTEM.4.6.0.0.station_files\\01\\010080']
In questa prima parte si andrà a ricalcare la procedura eseguita da CRUTEM per andare a creare la serie storica relativa alle anomalie di temperature nella zona della Groenlandia.
metadati = pd.DataFrame(columns=['ID','nome','paese','altezza','lat','lon'])
# creazione di un asse temporale che va da gennaio del 1850 al dicembre 2014
asset = pd.date_range('1850-01', '2015-01', freq='M')
nmonths=len(asset)
nyears=nmonths/12
#dati dalle stazioni.
dati = pd.DataFrame({'time': asset})
dati = dati.set_index(['time'])
# Come periodo di riferimento prendiamo il periodo tra il 1961 e il 1990
yref0 = 1961
yref1 = 1990
# caricamento delle osservazioni dalle stazioni
nst_tmp = nst
nodatacount=0
tooshortcount=0
outsidecount=0
from tqdm import tqdm
for si in tqdm(range(0,nst_tmp)):
filein=flist[si]
with open(filein) as f: data = f.readlines()
skipi = data.index("Obs:\n")+1
if (len(data)-skipi < 2):
nodatacount+=1
continue
yr = np.genfromtxt(filein, skip_header=skipi, delimiter=None, usecols=0, dtype='i4')
nyr = len(yr)
# primo filtro (se < di 30 anni scarto il file)
if (nyr < 30):
tooshortcount+=1
continue
# secondo filtro (country), prendiamo solo le stazioni relative alla Groenlandia
stcountry_line = [line for line in data if "Country=" in line]
stcountry = str(stcountry_line).split("=",1)[1]
stcountry = stcountry[:-4].strip()
if (stcountry != 'GREENLAND'):
continue
#terzo filtro (se presenti valori null)
norm_line = [line for line in data if "Normals=" in line]
norm = str(norm_line).split("=",1)[1]
norm = norm[:-4].strip()
if (norm.find("-99.0") > -1):
outsidecount+=1
continue
stname_line = [line for line in data if "Name=" in line]
stname = str(stname_line).split("=",1)[1]
stname = stname[:-4].strip()
country_line = [line for line in data if "Country=" in line]
country = str(country_line).split("=",1)[1]
country = country[:-4].strip()
elev_line = [line for line in data if "Height=" in line]
elev = str(elev_line).split("=",1)[1]
elev = float(elev[:-4].strip())
lat_line = [line for line in data if "Lat=" in line]
lat = str(lat_line).split("=",1)[1]
lat = float(lat[:-4].strip())
lon_line = [line for line in data if "Long=" in line]
lon = str(lon_line).split("=",1)[1]
lon = float(lon[:-4].strip())
metadati = metadati.append({
"ID": filein[-6:],
"nome": stname,
"paese": country,
"altezza": elev,
"lat": lat,
"lon": lon,
}, ignore_index=True)
ti0 = datetime.strptime(str(yr[0]), "%Y")
ti1 = datetime.strptime(str(yr[nyr-1]+1), "%Y")
stime = pd.date_range(ti0, ti1, freq='M')
loc_data_2d = np.genfromtxt(filein, skip_header=skipi, filling_values='-99.0', delimiter=None, dtype='float')[:,1:13]
loc_data_1d = loc_data_2d.flatten()
xdf = pd.DataFrame({'time': stime, filein[-6:]: np.asarray(loc_data_1d)})
xdf = xdf.set_index('time')
dati = pd.merge(dati, xdf, on='time', how='left')
10%|█ | 1037/10295 [00:11<01:43, 89.83it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
10%|█ | 1047/10295 [00:11<02:12, 69.62it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
10%|█ | 1055/10295 [00:11<02:46, 55.66it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
10%|█ | 1062/10295 [00:11<03:12, 47.95it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
10%|█ | 1068/10295 [00:12<03:44, 41.15it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
10%|█ | 1075/10295 [00:12<03:19, 46.19it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
11%|█ | 1081/10295 [00:12<03:45, 40.84it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
11%|█ | 1086/10295 [00:12<03:41, 41.62it/s]C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
16%|█▌ | 1657/10295 [00:18<01:28, 97.82it/s] C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\1659791162.py:67: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
metadati = metadati.append({
100%|██████████| 10295/10295 [00:50<00:00, 204.78it/s]
# stazioni prese in considerazione
metadati
| ID | nome | paese | altezza | lat | lon | |
|---|---|---|---|---|---|---|
| 0 | 042020 | PITUFFIK/THULE USAF | GREENLAND | 70.0 | 76.5 | 68.8 |
| 1 | 042100 | UPERNAVIK | GREENLAND | 63.0 | 72.8 | 56.2 |
| 2 | 042110 | UPERNAVIK | GREENLAND | 63.0 | 72.8 | 56.2 |
| 3 | 042180 | QEQERTARSUAQ | GREENLAND | 24.0 | 69.2 | 53.5 |
| 4 | 042200 | EGEDESMINDE | GREENLAND | 47.0 | 68.7 | 52.8 |
| 5 | 042210 | ILULISSAT/JAKOBSHAVN | GREENLAND | 31.0 | 69.2 | 51.1 |
| 6 | 042300 | SISIMIUT | GREENLAND | 12.0 | 66.9 | 53.7 |
| 7 | 042310 | KANGERLUSSUAQ (SDR. | GREENLAND | 53.0 | 67.0 | 50.7 |
| 8 | 042500 | GODTHAB/NUUK | GREENLAND | 50.0 | 64.2 | 51.8 |
| 9 | 042600 | PAAMIUT (FREDERIKSH | GREENLAND | 15.0 | 62.0 | 49.7 |
| 10 | 042700 | MITTARFIK/NARSARSUAQ | GREENLAND | 31.0 | 61.2 | 45.4 |
| 11 | 042720 | JULIANEHAB/QAQORTOQ | GREENLAND | 32.0 | 60.7 | 46.1 |
| 12 | 043100 | NORD | GREENLAND | 35.0 | 81.5 | 16.8 |
| 13 | 043200 | DANMARKSHAVN | GREENLAND | 11.0 | 76.8 | 18.8 |
| 14 | 043300 | DANEBORG | GREENLAND | 44.0 | 74.3 | 20.2 |
| 15 | 043390 | SCORESBYSUND | GREENLAND | 65.0 | 70.5 | 22.0 |
| 16 | 043400 | KAP TOBIN | GREENLAND | 41.0 | 70.5 | 22.0 |
| 17 | 043500 | APUTITEEQ | GREENLAND | 20.0 | 67.8 | 32.3 |
| 18 | 043600 | AMMASALIK | GREENLAND | 50.0 | 65.6 | 37.6 |
| 19 | 043900 | PRINS CHRISTIANS S. | GREENLAND | 76.0 | 60.0 | 43.1 |
| 20 | 099999 | SW COMBINED SERIES | GREENLAND | -999.0 | -99.9 | -999.9 |
Prendendo in considerazione solamente la Groenlandia sono state trovate un totale di 20 stazioni, notiamo che l'ultima stazione con ID 099999 non è ben definita e quindi la eliminiamo.
metadati = metadati.drop(labels=[20], axis=0)
dati = dati.drop(labels=['099999'], axis=1)
Conversione dei missing value in NaN
dati[dati == -99.0] = np.nan
metadati = metadati.set_index(["ID"])
# correzione della longitudine
metadati.lon = -metadati.lon
Verifichiamo che tutte le stazioni prese in considerazione siano posizionate nella Groenlandia.
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.stock_img()
ax.scatter([metadati.lon],[metadati.lat],color='r',marker='o',s=1.0)
plt.show()
Verifichiamo stazione per stazione le temperature assolute.
fig = plt.figure(figsize=(13,19))
no = 20
subplots = (no,2)
n_panels = subplots[0] * subplots[1]
proj = ccrs.PlateCarree()
ssite = metadati.index
for fi, f in enumerate(ssite):
ax = fig.add_subplot(subplots[0], subplots[1], (fi*2)+1, projection=proj)
ax.set_title(' - '.join([metadati.nome[ssite[fi]],metadati.paese[ssite[fi]]]))
ax.stock_img()
ax.coastlines()
plt.plot(metadati.lon[ssite[fi]], metadati.lat[ssite[fi]],
color='red', marker='o', markersize=5,
transform=ccrs.Geodetic())
tser = fig.add_subplot(subplots[0], subplots[1], (fi+1)*2)
plt.plot(dati[ssite[fi]])
fig.tight_layout()
plt.show()
Dal plot vediamo tutte le stazioni che hanno registrato dati relativi alle temperature in Groenlandia. Ciascuna stazione ha cominciato a registrare i dati in anni diversi, questo significa che non si può avere una copertura totale delle temperature, inoltre quasi tutte le stazioni presentano dei missing value, consideriamo che stiamo studiando un fenomeno fisico si è deciso di non andare a sostituire i missing value, questo perché la temperatura varia per via di diversi fattori in gioco e quindi sarebbe necessario uno studio più approfondito in tal senso.
Modifiche effettuate:
metadati.reset_index(inplace = True)
metadati = metadati.drop(labels=[1], axis=0)
dati = dati.drop(labels=['042100'], axis=1)
metadati = metadati.drop(labels=[3], axis=0)
dati = dati.drop(labels=['042180'], axis=1)
metadati = metadati.set_index(["ID"])
x1=dati.index.get_loc('1961-01-31')
x2=dati.index.get_loc('1990-12-31')
# calcolo delle normals, si è fatta la media di tutti i gennaio dal 1960 al 1990 per la stazione 042020 e così via
# per tutti i mesi e per tutte le stazioni
data_normals = dati[x1:x2].groupby(dati[x1:x2].index.month).mean()
print(data_normals)
042020 042110 042200 042210 042300 042310 \
time
1 -23.580000 -17.956667 -13.363333 -13.373333 -12.733333 -18.889286
2 -25.193333 -20.403333 -15.576667 -14.876667 -13.850000 -19.000000
3 -23.170000 -20.373333 -16.246667 -15.150000 -13.990000 -17.660714
4 -16.103333 -13.263333 -9.613333 -8.093333 -7.120000 -7.757143
5 -4.496667 -3.900000 -1.820000 0.263333 -0.256667 2.675000
6 2.226667 1.743333 2.683333 4.563333 3.490000 8.539286
7 5.146667 5.313333 5.723333 7.570000 6.213333 10.789286
8 4.226667 5.030000 5.320000 6.426667 6.044828 8.625000
9 -2.180000 0.866667 2.353333 2.613333 3.250000 3.278571
10 -10.586667 -4.246667 -2.216667 -3.736667 -1.780000 -5.539286
11 -17.426667 -8.866667 -5.903333 -7.973333 -5.826667 -11.853571
12 -21.741379 -14.055172 -9.696552 -11.172414 -9.875862 -16.262963
042500 042600 042700 042720 043100 043200 \
time
1 -7.430000 -6.586667 -6.776667 -5.493333 -30.292000 -23.113333
2 -7.786667 -6.406667 -6.123333 -5.036667 -30.512500 -24.270000
3 -7.973333 -5.993333 -5.113333 -4.406667 -29.928000 -23.353333
4 -3.836667 -2.330000 -0.133333 -0.586667 -23.157692 -17.160000
5 0.643333 1.426667 5.151724 3.323333 -9.996000 -6.500000
6 3.920000 3.713333 8.337931 5.193333 -0.630769 0.796667
7 6.536667 5.553333 10.263333 7.173333 3.300000 3.717241
8 6.093333 5.310000 9.280000 7.220000 2.124000 2.303571
9 3.470000 3.573333 5.486667 4.983333 -8.368000 -4.327586
10 -0.683333 0.120000 0.363333 1.150000 -19.515385 -13.741379
11 -3.650000 -2.823333 -3.233333 -1.916667 -25.746154 -19.943333
12 -6.110345 -5.306897 -6.041379 -4.358621 -27.815385 -21.948276
043300 043390 043400 043500 043600 043900
time
1 -20.150000 -16.080000 -16.403846 -10.513043 -7.493333 -4.057692
2 -21.544000 -17.113333 -17.707692 -11.229167 -7.693333 -3.922222
3 -20.738462 -16.526667 -17.250000 -11.368000 -8.146667 -3.789655
4 -14.688889 -11.190000 -11.726923 -7.083333 -3.976667 -0.932143
5 -5.069231 -3.466667 -3.769231 -1.737500 0.673333 1.803448
6 1.160000 1.123333 0.576000 1.008333 4.180000 4.373333
7 4.004167 3.348276 2.734615 2.208696 6.393333 6.526667
8 3.462500 3.514286 3.264000 2.356522 6.003333 6.500000
9 -2.226087 -0.372414 -0.400000 0.447826 2.993333 4.425000
10 -10.931818 -6.379310 -6.380769 -2.579167 -0.870000 1.339286
11 -17.118182 -12.213333 -12.192308 -7.200000 -4.753333 -1.260714
12 -18.895652 -14.844828 -15.128000 -9.728000 -7.286207 -3.134615
Ora dobbiamo calcolare le anomalie, e quindi sottrarre mese per mese a tutti i dati della serie temporale il corrispettivo mese delle climate normals stazione per stazione.
data_anom = dati
y0=0
for yi in range(0,int(nyears)):
data_anom.iloc[y0:y0+12,:] = dati.iloc[y0:y0+12,:]-np.array(data_normals.iloc[:,:])
y0+=12
fig = plt.figure(figsize=(13,19))
no = 20
subplots = (no,2)
n_panels = subplots[0] * subplots[1]
proj = ccrs.PlateCarree()
ssite = metadati.index
for fi, f in enumerate(ssite):
ax = fig.add_subplot(subplots[0], subplots[1], (fi*2)+1, projection=proj)
ax.set_title(' - '.join([metadati.nome[ssite[fi]],metadati.paese[ssite[fi]]]))
ax.stock_img()
ax.coastlines()
plt.plot(metadati.lon[ssite[fi]], metadati.lat[ssite[fi]],
color='red', marker='o', markersize=5,
transform=ccrs.Geodetic())
tser = fig.add_subplot(subplots[0], subplots[1], (fi+1)*2)
plt.plot(data_anom[ssite[fi]])
fig.tight_layout()
plt.show()
Dai grafici in alto possiamo notare le anomalie di temperatura per ciascuna stazione.
# si crea una griglia utilizzando una risoluzione di 10 gradi
resol = 10
nlon = np.int(360/resol)
nlat = np.int(180/resol)
grlons = np.empty([nlon+1],dtype='float')
grlats = np.empty([nlat+1],dtype='float')
grlons[0] = -180.
grlats[0] = -90.
for i in range(1,nlon+1):
grlons[i]=grlons[i-1]+resol
for i in range(1,nlat+1):
grlats[i]=grlats[i-1]+resol
C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\4271657429.py:3: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations nlon = np.int(360/resol) C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\4271657429.py:4: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information. Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations nlat = np.int(180/resol)
data_mo = np.empty([nmonths,nlat,nlon],dtype=float)
data_mo[:,:,:] = np.nan
data_yr = np.empty([int(nyears),nlat,nlon],dtype=float)
data_yr[:,:,:] = np.nan
# si va ad analizzare una alla volta in sequenza tutte le celle della griglia creata
for j in range(0,nlat):
for i in range(0,nlon):
# si verificano le coordinate delle singole stazioni andando a considerare solo quelle che sono all'interno della cella
dummy = metadati[metadati.lon >= grlons[i]]
dummy = dummy[dummy.lon < grlons[i+1]]
dummy = dummy[dummy.lat >= grlats[j]]
dummy = dummy[dummy.lat < grlats[j+1]]
# nel caso dentro la cella sia presente almeno una stazione si eseguono i passaggi successivi
if (len(dummy.index) > 0):
#calcolo della media delle serie temporali delle stazioni che sono all'interno della cella
data_mo[:,j,i] = np.array(data_anom[dummy.index].mean(axis=1)).flatten()
data_yr[:,j,i] = np.array(data_anom[dummy.index].mean(axis=1).resample("Y").mean()).flatten()
fig = plt.figure(figsize=(13,3))
subplots = (1,2)
n_panels = subplots[0] * subplots[1]
ax = fig.add_subplot(subplots[0], subplots[1], 1, projection=ccrs.PlateCarree())
ax.stock_img()
ax.coastlines()
ax.gridlines(xlocs=range(-180,181,resol), ylocs=range(-90,91,resol))
ax.scatter([dummy.lon],[dummy.lat],color='r',marker='o',s=1.2)
tser = fig.add_subplot(subplots[0], subplots[1], 2)
plt.plot(asset,np.array(data_anom[dummy.index]),color='grey',linewidth=0.1)
plt.plot(asset,data_mo[:,j,i],color='blue',linewidth=0.3)
plt.plot(pd.date_range('1850-01', '2015-01', freq='Y'),data_yr[:,j,i],color='red',linewidth=2)
fig.tight_layout()
plt.show()
Abbiamo ora le anomalie di temperature raggruppate per celle. Per ultimo si devono aggregare tutte le celle nello spazio, in modo da ottenere un'unica serie di anomalie che sia rappresentativa per l'intera Groenlandia.
fig = plt.figure(figsize=(13,13))
subplots = (2,1)
ax = fig.add_subplot(subplots[0], subplots[1], 1, projection=ccrs.PlateCarree())
ax.stock_img()
ax.coastlines()
ax.gridlines(xlocs=range(-180,181,resol), ylocs=range(-90,91,resol))
tser = fig.add_subplot(subplots[0], subplots[1], 2)
data_glob = np.empty([int(nyears)],dtype=float)
data_glob[:] = np.nan
for j in range(0,nlat):
for i in range(0,nlon):
if not(np.isnan(data_yr[:,j,i]).all()):
ax.plot( grlons[i:i+2].sum()/2, grlats[j:j+2].sum()/2,
color='red', marker='x', markersize=5,transform=ccrs.PlateCarree())
plt.plot(pd.date_range('1850-01', '2015-01', freq='Y'),data_yr[:,j,i],color='gray',linewidth=0.1)
for t in range(0,int(nyears)):
local = data_yr[t,:,:]
valid = np.isnan(data_yr[t,:,:])
data_glob[t] = np.nanmean(local[~valid])
plt.plot(pd.date_range('1850-01', '2015-01', freq='Y'),data_glob[:],color='red',linewidth=2)
plt.axhline(y=0, color='black', linestyle='-')
plt.title('Anomalie di temperatura in Groenlandia')
plt.ylabel('[°C]')
fig.tight_layout()
Anomalie di temperatura della Groenlandia.
df1 = pd.DataFrame(pd.date_range('1850', '2015', freq='Y'), columns = ['anno'])
df1['anno'] = pd.DatetimeIndex(df1['anno']).year
df1['anomGroe'] = data_glob
Per la creazione dell'AMO Index utilizziamo un modello relativo alle temperature superficiali in gradi C° dei mari di tutto il globo. Il periodo del modello va dal gennaio 1850 sino al dicembre 2014.
L'indice AMO viene calcolato prendendo la media globale annuale delle anomalie di temperature superficiale dei mari escludendo il Nord dell'Atlantico (0°-60°N, 75°W-7.5°W) sottraendola dalla media delle anomalie annuali calcolate esclusivamente sul Nord Atlantico.
modfile='./tos_Omon_MCM-UA-1-0_historical_r1i1p1f2_gn_185001-201412.nc'
Si comincia costruendo un modello che tiene conto dell'intero pianeta tranne il Nord dell'Altantico e poi un modello che tiene conto esclusivamente di quest'ultimo.
d1 = d2 = d3 = d4 = xr.open_dataset(modfile)
d1
<xarray.Dataset>
Dimensions: (bnds: 2, longitude: 192, latitude: 80, time: 1980)
Coordinates:
* bnds (bnds) float64 1.0 2.0
* longitude (longitude) float64 -0.9375 0.9375 2.812 ... 353.4 355.3 357.2
* latitude (latitude) float64 -88.63 -86.13 -83.88 ... 83.88 86.13 88.63
* time (time) object 1850-01-17 00:00:00 ... 2014-12-17 00:00:00
Data variables:
lon_bnds (longitude, bnds) float64 -1.875 0.0 0.0 ... 356.2 356.2 358.1
lat_bnds (latitude, bnds) float64 -90.0 -87.26 -87.26 ... 87.26 87.26 90.0
tos (time, latitude, longitude) float32 ...
time_bnds (time, bnds) object 1850-01-01 12:00:00 ... 2015-01-01 12:00:00
Attributes: (12/45)
source: Manabe Climate Model at U of Arizona ...
activity_id: CMIP
branch_method: standard
Conventions: CF-1.7 CMIP-6.2
creation_date: 2019-08-14T00:00:00Z
data_specs_version: 01.00.28
... ...
sub_experiment: none
sub_experiment_id: none
title: UArizona MCM-UA-1-0 historical
product: model-output
institution: Department of Geosciences, University of Arizona,...
license: CMIP6 model data produced by the U of Arizona is ...# conversione in datetime64
datetimeindex = d1.indexes['time'].to_datetimeindex()
datetimeindex = d2.indexes['time'].to_datetimeindex()
datetimeindex = d3.indexes['time'].to_datetimeindex()
datetimeindex = d4.indexes['time'].to_datetimeindex()
d1['time'] = datetimeindex
d2['time'] = datetimeindex
d3['time'] = datetimeindex
d4['time'] = datetimeindex
C:\Users\fadda\anaconda3\envs\time_series\lib\site-packages\xarray\coding\times.py:360: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self. sample = dates.ravel()[0] C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\4293979708.py:2: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar. This may lead to subtle errors in operations that depend on the length of time between dates. datetimeindex = d1.indexes['time'].to_datetimeindex() C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\4293979708.py:3: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar. This may lead to subtle errors in operations that depend on the length of time between dates. datetimeindex = d2.indexes['time'].to_datetimeindex() C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\4293979708.py:4: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar. This may lead to subtle errors in operations that depend on the length of time between dates. datetimeindex = d3.indexes['time'].to_datetimeindex() C:\Users\fadda\AppData\Local\Temp\ipykernel_15492\4293979708.py:5: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar. This may lead to subtle errors in operations that depend on the length of time between dates. datetimeindex = d4.indexes['time'].to_datetimeindex()
Costruzione del modello relativo alla parte Nord Atlantica
modelloAmo = d1.where((d1.latitude > 0) & (d1.latitude < 60) & (d1.longitude > 285) & (d1.longitude < 352.5), drop=True)
p0 = modelloAmo.tos.mean('time').plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet',extend='both')
p0.axes.coastlines()
p0.axes.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x2124d8bbd30>
Costruzione del modello escludendo l'Atlantico del Nord
modello2 = d2.where((d2.latitude < 90) & (d2.longitude < 352.5) & (d2.longitude <= 285) , drop=False)
p0 = modello2.tos.mean('time').plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet',extend='both')
p0.axes.coastlines()
p0.axes.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x212513d9e50>
modello3 = d3.where((d1.latitude < 0) & (d3.longitude > 290) & (d3.longitude < 360), drop=False)
p0 = modello3.tos.mean('time').plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet',extend='both')
p0.axes.coastlines()
p0.axes.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x2124f851130>
modello4 = d4.where((d1.latitude >= 60), drop=False)
p0 = modello4.tos.mean('time').plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet',extend='both')
p0.axes.coastlines()
p0.axes.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x2124d917ee0>
Si effettua un merge del modello2, modello3 e modello4, in modo da ricostruire il modello di interesse per noi.
modelloNoAtl = xr.merge([modello2,modello3,modello4])
p0 = modelloNoAtl.tos.mean('time').plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet',extend='both')
p0.axes.coastlines()
p0.axes.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x2124fdfe550>
Dopo aver costruito i due modelli andiamo a calcolare le anomalie di temperatura:
# Calcoliamo le normals prendendo il periodo di riferimento, raggruppando per mese e andando a effettuare la media
normalNoAtl = modelloNoAtl.tos.sel(time=slice('1850-01', '1890-12-31')).groupby('time.month').mean()
normalAmo = modelloAmo.tos.sel(time=slice('1850-01', '1890-12-31')).groupby('time.month').mean()
# Le anomalie vengono calcolate raggruppando per mese le temperature assolute e poi sottraendo le normals
anomalieNoAtl = modelloNoAtl.tos.groupby('time.month') - normalNoAtl
anomalieAmo = modelloAmo.tos.groupby('time.month') - normalAmo
# modifica dei pesi
# Visto che le aree delle celle diminuiscono con la latitudine si vanno a creare dei pesi proporzionali alle celle
# in modo da poter creare una media pesata
pesiNoAtl = np.cos(np.deg2rad(modelloNoAtl.latitude))
pesiAmo = np.cos(np.deg2rad(modelloAmo.latitude))
# Il modello completo conterrà le anomalie di temperature dell'Atlantico del Nord e del globo intero escludendo il Nord Atlantico.
# Da una risoluzione mensile passiamo ad una risoluzione annuale, aggregando poi anche nello spazio.
modelloCompleto = xr.Dataset()
modelloCompleto['tNoAtl'] = anomalieNoAtl.weighted(pesiNoAtl).mean(('longitude', 'latitude'), keep_attrs=True).resample(time='Y').mean()
modelloCompleto['tAmo'] = anomalieAmo.weighted(pesiAmo).mean(('longitude', 'latitude'), keep_attrs=True).resample(time='Y').mean()
Seguendo la definizione dell'indice AMO, si va quindi a sottrarre alle anomalie dell'Atlantico quelle relative a tutto il globo escludendo il Nord Atlantico
totale = modelloCompleto.tAmo - modelloCompleto.tNoAtl
totale
<xarray.DataArray (time: 165)>
array([-0.07263093, -0.22366853, -0.28658982, -0.0055107 , -0.22068568,
-0.15473103, -0.26443282, -0.14798923, -0.21817076, -0.221023 ,
-0.20454415, -0.11603901, 0.06204571, 0.15235264, 0.01780326,
-0.04170245, -0.18067771, -0.08880673, 0.13044599, -0.21868026,
0.06800493, 0.19949009, 0.16899327, 0.56415666, 0.21012964,
0.04514276, -0.09230632, 0.11883762, 0.14837877, 0.16980886,
-0.11246507, 0.18595026, 0.18903891, 0.14902244, 0.02561967,
-0.11485682, -0.05632997, 0.18056069, -0.15458829, 0.17137712,
0.23927074, 0.04803847, 0.03627353, 0.2803443 , 0.11586739,
-0.01859957, -0.08989961, -0.32751925, 0.00724594, 0.00996242,
0.17878007, 0.10384548, 0.10215164, 0.0776101 , 0.01139102,
0.01677712, 0.06992126, 0.15792389, 0.11136409, 0.11155269,
0.040364 , -0.1133489 , -0.3970292 , -0.23458152, -0.22976045,
0.09084673, -0.10571517, -0.00772636, 0.0859891 , -0.00574517,
0.04027676, 0.08953073, -0.11436532, -0.20162388, -0.24563517,
0.04430177, 0.14587949, -0.08706278, 0.0955557 , 0.14624774,
-0.11901757, 0.06409697, -0.01800325, 0.04904631, 0.01416419,
0.18627258, 0.0519064 , -0.05746475, -0.04584225, 0.14548027,
0.02536201, 0.04335942, 0.16644177, -0.19813158, -0.07018882,
-0.21742252, -0.07218918, -0.16799216, -0.32103125, -0.12199013,
0.02659784, -0.22993997, 0.13123566, -0.14043248, -0.05923033,
0.07452597, 0.01401333, 0.10343263, 0.00709666, -0.1789199 ,
-0.05340617, -0.34674803, -0.3652285 , -0.23773795, -0.11978361,
-0.44787504, -0.32234429, -0.3513176 , -0.12119485, 0.04117268,
-0.0026219 , -0.21733221, -0.23179569, -0.06486765, -0.04235908,
0.15465302, 0.06245346, -0.10211512, -0.18674023, -0.20051849,
-0.33691783, -0.19193922, -0.26459378, -0.3018955 , -0.30659096,
-0.20916768, 0.01497165, -0.00269398, -0.0480597 , 0.16328518,
0.21509416, 0.01457634, 0.02589139, -0.03650928, 0.06242263,
0.26190217, 0.12469045, 0.16163021, 0.13165462, 0.12400276,
0.15535545, -0.08205654, -0.1053201 , -0.19202144, -0.01987248,
0.20346015, 0.07531302, 0.08818656, 0.13871626, -0.10674908,
0.18720088, -0.08335551, -0.31027143, -0.34111918, -0.34566067])
Coordinates:
* time (time) datetime64[ns] 1850-12-31 1851-12-31 ... 2014-12-31Costruiamo un dataframe che conterrà anno per anno il valore dell'AMO Index.
df2 = pd.DataFrame(pd.date_range('1850', '2015', freq='Y'), columns = ['anno'])
df2['anno'] = pd.DatetimeIndex(df2['anno']).year
df2['anomAmo'] = totale
df2 = df2.set_index('anno')
plt.figure(figsize=(13,8))
x = df2.anomAmo
y = np.linspace(1860, 2020, 165)
plt.axhline(y=0, color='black', linestyle='-')
plt.title('Oscillazione multidecennale atlantica')
plt.ylabel('AMO Index')
plt.fill_between(y, x, where=x>0, color='orange')
plt.fill_between(y, x, where=x<0, color='blue')
plt.show()
Nel grafico sopra si può osservare l'andamento di questo indice, e quindi le sue oscillazioni sia negative che positive dal 1860 al 2014.
In basso vediamo come dovrebbe risultare l'AMO Index.

df2.reset_index(inplace = True)
finale = df1.merge(df2, how='left')
finale = finale.set_index('anno')
plt.figure(figsize=(12,5))
ax3 = finale.anomAmo.rolling(window =15).mean().plot(color='green', grid=True, label='AMO Index')
ax4 = finale.anomGroe.rolling(window =15).mean().plot(color='red', grid=True, secondary_y=True, label='Anomalie Groenlandia')
plt.title('Confronto tra le anomalie di temperatura in Groenlandia e AMO Index')
plt.xlabel('Anno')
ax3.set_ylabel('AMO Index', color='green')
ax4.set_ylabel('Anomalie Groenlandia [°C]', color='red')
h1, l1 = ax3.get_legend_handles_labels()
h2, l2 = ax4.get_legend_handles_labels()
plt.show()
Nel grafico sopra sono state plottate due serie storiche: in verde l'AMO Index e in rosso le anomalie di temperatura nella Groenlandia, per una migliore visualizzazione si è usata una media mobile per entrambe le serie.
Osservando il grafico si può osservare come non ci sia una netta tendenza tra le due, ma è vero che nel periodo tra il 1925 - 1940 l'incremento dell'AMO Index si riflette anche in un aumento delle anomalie di temperatura della Groenlandia, e allo stesso modo si nota un decremento sia delle anomalie che dell'AMO Index nei decenni successivi. Si nota infine, un andamento divergente negli anni 2000 tra l'AMO Index e le anomalie di temperatura, questo potrebbe essere influenzato dall'aumento delle anomalie di temperature globali che si sono registrate a partire del 21° secolo.
Per concludere si può quindi affermare che da una semplice visualizzazione grafica l'AMO Index ha una certa influenza sulle temperature della Groenlandia, nell'ultimo decennio però l'andamento divergente che si nota dal grafico potrebbe essere indicativo di un cambio di passo rispetto all'influenza che l'AMO esercita.